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4.  SCIENTIFIC  PROGRESS  AND  ACCOMPLISHMENT.  This  project  is  about  the  numerical 
and  theoretical  study  of  granular  flows.  Significant  advances  have  been  achieved  in  relation  to  numerical 
modeling.  People  who  have  collaborated  to  the  project  during  this  period  of  funding  include  Profs.  D.G. 
Schaeffer  (Duke  University)  and  M.  Shearer  (North  Carolina  State  University),  Dr.  J.V.  Matthews  (Duke 
University)  and  K.A  Coffey  (Ph.D.  student,  Gremaud  adviser).  We  also  work  in  consultation  with  Engineers 
at  Jenike  &  Johanson,  Inc.,  Westford,  MA.,  and  in  particular  T.A.  Royal,  Vice  President,  to  develop  efficient 
and  robust  solvers  for  various  types  of  problems  related  to  granular  flows. 

One  of  the  first  and  main  problems  one  encounters  when  dealing  with 
granular  materials  is  to  obtain  a  proper  mathematical  model.  Several  mod¬ 
els  exist,  that  appear  to  capture  well  some,  but  in  general  not  all,  aspects 
of  the  phenomena  under  consideration.  Our  approach  is  twofold.  First,  we 
aim  at  obtaining,  through  careful  Analysis  and  Numerical  Analysis  of  ex¬ 
isting  models,  efficient  and  reliable  computational  tools  that  will  allow  for 
meaningful  investigations  the  models  themselves,  through  comparisons  with 
experiments  for  instance.  Second,  and  at  a  less  fundamental  level,  having 
such  computational  tools  represent  a  breakthrough  in  this  field,  and  open 
new  opportunities  in  terms  of  numerical  simulations  in  industrial  settings. 

Three  main  fields  of  applications  have  been  considered:  hopper  flows, 
powder  consolidation  and  sound  propagation  in  granular  materials.  Our 
main  contributions  during  the  last  year  have  been 

•  identification  of  secondary  circulation  phenomena  and  resonance  in  three  dimensional  hoppers;  this 
observation  is  astonishing  and  could  potentially  revolutionize  silo  design; 
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Fig.  4.2.  Stream  function  showing  secondary  flow;  top:  tilted  hopper  (m  =  1),  bottom:  ’’elliptical”  hopper  (m  =  2).  Angle 
of  wall  friction  =  15°  (n  =  tan  15°),  half  opening  angle  6W  =  30°.  By  symmetry,  only  half  of  the  hopper  is  represented. 


•  further  progress  in  the  implementation  (Objective-C)  of  our  model  for  powder  consolidation  to  ax- 
isymmetric  domains  of  varying  cross  section  in  an  industrial  setting, 

•  numerical  analysis  of  sound  propagation  phenomena  in  granular  materials;  our  analysis, the  first  one 
done  using  higher  order  methods  (third  order  Discontinuous  Galerkin  method)  solves  a  controversial 
question  which  had  been  incorrectly  answered  in  the  literature  as  a  result  of  substandard  numerics. 

We  now  proceed  to  describe  in  more  details  the  first  and  third  points. 

4.1.  Secondary  circulation  in  hopper  flows.  We  have  studied  3D  steady  flow  in  a  pyramidal  hopper 

^  =  {(r,6>,0)  :0  <<9<C(0)}. 

This  work  prepares  the  way  to  study  flow  in  fully  three-dimensional  geometries,  such  as  the  transition  hopper 
illustrated  in  Figure  4.1  or  a  silo  with  multiple  outlets.  We  have  begun  the  study  of  flows  using  a  perturbation- 
theory  analysis.  For  conical  hoppers  (i.e.,  C(0)  =  6W,  a  constant),  classical  similarity  solutions  can  be 
constructed.  Those  so-called  Jenike  solutions  are  such  that  only  the  radial  component  of  the  velocity  is 
nonzero,  all  variables  are  independent  of  0,  and  the  r-dependence  of  the  solution  has  the  similarity  form 

vr  =  —r~2  v(0),  T  =  r  T(0). 

Suppose  the  function  C  specifying  the  boundary  of  Q  has  the  expansion 

C(4>)  =  6W  +  e  cos (m<j))  +  0(e2); 

for  example,  a  slightly  tilted  (circular)  cone  admits  such  an  expansion  with  m  —  1,  where  e  measures  the 
angle  of  tilt;  likewise  for  a  (vertical)  pyramidal  hopper  having  a  slightly  elliptical  cross  section,  with  m  =  2. 
An  expansion  of  the  solution 


v  =  u(0)  +  eu(1)  +  0(e2),  T  =  T^0)  +  eT(1)  +  0(e2), 

can  be  sought,  where  have  the  form  of  the  Jenike  solution.  A  numerical  code  based  on  a  differential- 

algebraic  formulation  has  been  written  to  solve  the  boundary  value  problem  one  obtains  for  T^.  This 
analysis  has  brought  to  the  fore  some  extremely  interesting  phenomena.  In  particular  all  three  components 
of  v need  to  be  nonzero:  in  other  words,  secondary  circulation  appears  when  axial  symmetry  is  broken! 
This  circulatory  flow  is  illustrated  in  Figure  4.2  (left)  for  a  tilted  hopper  m  =  1  and  Figure  4.2  (right)  for  an 
elliptical  hopper  m  =  2.  In  other  words,  the  grains  do  not  move  along  radial  lines,  as  is  the  case  for  Jenike’ s 
solutions  but  follow  more  complicated  and  fully  three  dimensional  trajectories.  Very  singular  behavior  occurs 
when  the  parameters  in  the  problem  are  varied.  For  example,  consider  Figure  4.3  (left),  which  shows,  for 
m  «  1  and  other  parameters  fixed,  the  circulation  velocity  as  a  function  of  the  opening  angle  of  the 

hopper:  suffers  a  “l/x-blow-up”  as  0W  passes  through  a  critical  value. 

While  this  blow  up  could  not  have  been  anticipated,  a  posteriori ,  its  cause  is  easily  understood:  the  two- 
point  boundary  problem  for  has  inhomogeneous  boundary  conditions,  but  at  the  singular  point  in 
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as  a  function  of  the  wall  friction  pw  (internal  friction  S  =  30°,  half  opening  angle  6W  =  30°).  Right:  critical  values  leading  to 
sign  changes  of  the  circulation  (internal  friction  S  =  30°,  m  =  1). 


Figure  4.3  (left),  the  problem  with  homogeneous  boundary  conditions  has  a  nonzero  solution.  This  behavior 
is  analogous  to  that  of  a  forced  linear  harmonic  oscillator 


x  +  (JqX  =  A  sin  cat  : 


the  steady-state  oscillations  have  amplitude  A/(uJq  —  cj2),  which  diverges  as  u  passes  through  the  natural 
frequency  cjo  of  the  oscillator.  Unlike  the  above  problem,  in  our  model  the  inhomogeneity  is  in  the  equation 
rather  than  a  boundary  condition.  Based  on  this  analogy,  we  describe  this  divergence  as  a  kind  of  resonance. 

As  0W  crosses  the  critical  value  in  Figure  4.3  (left),  the  direction  of  circulation  changes.  It  turns  out 
that  there  are  also  additional  parameter  values  for  which  the  sign  of  the  circulation  changes  smoothly ,  passing 
continuously  through  zero.  Curves  of  0W,  fiw  along  which  the  circulation  changes  sign  by  either  mechanism  are 
shown  in  Figure  4.3  (right).  Those  are  fundamentally  new  observations.  They  lead  to  fascinating  questions 
discussed  in  Section  4.3.1. 

4.2.  Wave  propagation  in  granular  materials.  The  time  dependent  versions  of  most  plasticity  mod¬ 
els  for  granular  flows  are  hopelessly  ill-posed.  We  have  investigated  in  some  simple  cases  how  equilibrium  is 
reached.  More  precisely,  considering  the  granular  equivalent  of  the  waterhammer  problem,  see  Figure  4.4,  it 
has  been  argued  by  some  authors1  that  the  classical  “Janssen  stress  profile”  for  a  granular  material  at  rest  is 
not  the  asymptotic  solution  that  is  reached  at  equilibrium.  Considering  the  experiment  depicted  in  Figure  4.4, 
the  two  main  unknowns  are  u  and  u,  the  strain  and  the  velocity.  If  a  stands  for  the  average  horizontal  stress, 
the  equations  in  Lagrangian  coordinates  take  the  form 

(4.1)  dtu-dxv  =  0, 

(4.2)  dtv  +  dxd  G  —k  sgn(u)  cr, 

where  k  is  a  material  parameter  depending  on  both  the  internal  friction  of  the  material  and  the  frictional 
effects  between  the  wall  and  the  grains.  The  stress  a  is  taken  as  a  function  of  the  strain  u. 

The  propagation  of  acoustic  waves  in  dry  granular  materials  and  the  corresponding  influence  of  wall 
friction  can  then  be  studied  in  the  above  framework.  The  way  an  asymptotic  profile  is  reached  can  be 
numerically  determined  and  is  illustrated  in  Figure  4.5  which  shows  the  evolution  of  a  pressure  step  toward  a 


xSee  for  instance  T.  Boutreux,  E.  Raphael  and  P.G.  de  Gennes,  Propagation  of  a  pressure  step  in  a  granular  material:  the 
role  of  wall  friction ,  Phys.  Rev.  E,  55  (1997),  pp. 5759-5773. 
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Fig.  4.4.  “Waterhammer”  experiment:  a  semi  infinite  cylinder  contains  a  compressible  granular  material;  a  pressure  step  is 
applied  at  the  end. 


Fig.  4.5.  Propagation  of  a  pressure  wave  in  a  granular  material  with  a  quadratic  stress-srain  relationship  (m  =  2)  and 
convergence  toward  a  Janssen’s  like  asymptotic  state;  solutions  are  plotted  at  time  t  =  1,  2,  3,  5,  10  and  200.  Left:  stress,  right: 
velocity 


Janssen-like  state.  A  pressure  front,  followed  by  a  stopping  wave,  propagate  through  the  material.  When  the 
latter  catches  up  with  the  former,  all  motion  stops. 

The  above  system  is  a  hyperbolic  system  of  nonlinear  conservation  laws.  Using  a  high  order  (3)  Dis¬ 
continuous  Galerkin  method  together  with  a  careful  implementation  of  the  right-hand  term  corresponding  to 
a  graph,  we  have  been  to  show  that  the  solutions  to  the  above  system  do  indeed  convergence  to  Janssen’s 
solution.  We  explain  the  discrepancy  between  our  results  and  some  previously  published  results  by  other 
authors  by  the  excess  of  numerical  diffusion  introduced  by  their  methods. 

4.3.  Ongoing  and  future  work. 

4.3.1.  Hopper  flows.  The  main  progress  on  hopper  hows  during  this  period  of  funding  has  been  the 
discovery  of  secondary  circulations  in  most  pyramidal  hoppers.  We  do  not  argue  that  secondary  circulation 
actually  takes  in  lab  experiments.  We  claim,  however,  that  the  current  models  predict  them,  a  fact  that  had 
escaped  everybody  so  far,  mainly  because  a  full  numerical  study  would  require  the  numerical  resolution  of 
a  nonlinear  elliptic  DEA  (Differential  Algebraic  Equations)  system  of  9  unknowns.  We  are  in  close  contact 
with  Prof.  R.H.  Behringer  in  the  Physics  Department  at  Duke  who  is  interested  setting  up  experiments  that 
will  bring  to  the  fore  secondary  circulation.  In  case  circulation  phenomena  are  indeed  observed,  that  will 
add  a  fundamental  aspect  to  the  present  held.  Alternatively,  the  absence  of  such  phenomena  would  call  for  a 
complete  and  thorough  revision  of  the  whole  held. 

Should  the  experiments  support  the  theory,  we  plan  to  discretize  the  full  three  dimensional  problem  in 
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the  next  period  of  funding. 

4.3.2.  Powder  consolidation.  Our  previous  efforts  toward  the  implementation  of  a  reliable  and  efficient 
powder  consolidation  code  have  been  expanded.  Having  obtained  a  satisfactory  model,  we  have  started 
its  implementation  in  Objective-C  to  ease  the  inclusion  of  the  corresponding  code  into  existing  dedicated 
engineering  packages.  We  expect  to  have  a  first  version  ready  by  July  2002,  when  our  industrial  contact  Tony 
Royal  will  visit  us  for  10  days  this  summer. 

4.3.3.  Wave  propagation  in  granular  material.  The  problem  of  wave  propagation  as  set  in  Section 
4.2  is  interesting  in  several  respects. 

For  instance,  most  existing  models  assume  a  stress-strain  relationship  of  the  type  a(u)  =  —ao  sgn (u)  \u\m, 
<Jo  >  0.  However,  theoretical  and  experimental  studies  differ  as  to  what  the  value  of  m  is,  namely  3/2  and 
2,  respectively.  The  present  problem  is  an  ideal  testbed  to  settle  this  question.  Again,  our  colleague  Bob 
Behringer  has  expressed  strong  interest  in  directing  a  series  of  experiments  which, by  comparision  with  our 
numerical  results  will  allow  the  determination  of  m. 

Mathematically,  the  problem  is  also  very  interesting  and  open.  Closely  related  systems  with  damping 
have  recently  received  a  quite  a  lot  of  attention.  In  all  those  contributions  the  damping  terms  are  proportional 
to  the  velocity  instead  being  proportional  to  the  stress  a  as  is  the  case  for  dry  friction.  The  present  case  is 
much  tougher  to  analyze,  as  no  energy  estimate  can  be  easily  derived.  We  plan  on  investigating  As  part  of 
this  project  conditions  of  existence  of  weak  solutions  to  problems  of  the  present  type  will  be  sought. 

As  a  generalization  on  the  numerical  side,  we  will  also  develop  efficient  and  reliable  numerical  methods 
for  the  resolution  of  hyperbolic  systems  with  either  stiff,  discontinuous  or  set- valued  right  hand  side  will  be 
constructed,  studied  and  implemented. 

4.3.4.  PDAEs.  PDAs  are  ubiquitous  in  plasticity  theory.  One  of  our  main  long  term  goals  is  the 
development  of  numerical  methods  to  efficiently  handle  this  type  of  mathematical  structure.  Several  methods 
will  be  tested  including  the  use  ”  streamfunctions”  leading  to  systems  of  Hamilton- Jacobi  equations,  adaptation 
of  DAE  techniques  such  as  constraint  stabilization  methods  to  PDAEs,  use  of  variational  formulations  such 
as  constraint  least  gradients. 

4.3.5.  Frictional  chutes.  The  design  of  optimal  chutes  is  important  in  countless  applications.  We  have 
started  the  study  of  finding  optimal  three  dimensional  chutes  when  the  particles  are  subject  to  Coulomb’s  dry 
friction.  In  first  approximation,  the  chute  is  considered  as  a  wire  and  the  particles  as  beads  allowed  to  slide 
along  that  wire.  For  given  location  and  direction  of  the  starting  point  and  outlet  of  the  chute,  the  problem 
consists  in  finding  the  optimal  geometry  linking  the  two.  Optimal  is  does  not  refer  to  minimizing  time,  as  in 
the  classical  brachistochrone  problem,  but  rather  to  minimizing  the  work  done  by  friction.  A  typical  industrial 
application  is  the  handling  of  brittle  material  such  as  pasta  for  instance. 

4.3.6.  Granular  heap.  When  a  dry  granular  material  is  poured  on  a  surface,  the  slope  of  the  resulting 
pile  cannot  exceed  some  critical  value:  the  angle  of  repose.  The  pile  takes  a  conical  shape,  such  as  for  instance 
that  of  sand  in  the  lower  part  of  an  hour  glass.  In  this  project,  we  study  how  to  determine  the  shape  of  such 
piles  when  obstacles  are  present,  see  Figure  4.6,  right. 

A  deceptively  simple  problem  is  the  determination  of  the  volume  of  material  inside  a  bin  of  possibly 
complicated  geometry  including  inserts,  internal  walls,  etc...  In  practice,  this  is  done  by  having  one  or  several 
gauges  measuring  the  height  of  the  material  at  given  points.  Without  some  knowledge  of  the  geometry  of  the 
top  free  boundary,  such  methods  are  inaccurate.  This  problem  has  a  long  history,  largely  due  to  its  connection 
to  soil  mechanics  in  civil  and  military  engineering.  Coulomb  was  the  first  to  relate  the  angle  of  slip  to  the 
friction  properties  of  the  material.  Unlike  previous  work,  he  did  not  assume  a  priori  values  of  that  angle. 

We  are  proposed  a  new  Eikonal  formulation  of  those  problems.  Problems  with  cylindrical  obstacles  (of 
arbitrary  cross  section)  correspond  to  traveltime  problems  with  obstacles.  When  the  obstacle  is  not  cylindrical, 
the  formulation  is  more  involved,  but  still  Eikonal.  We  have  supposed,  studied  and  implemented  numerical 
methods  based  on  modification  of  the  fast  marching  method,  that  preserve  both  the  sweeping  (one  sweep 
instead  of  many  iterations!)  character  of  the  method  and  preserve  the  order  of  accuracy  (up  to  and  included 
second  order),  even  in  the  presence  of  obstacles  (domains  of  infinite  slowness). 
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Heap  with  6-40  and  source  at  (0,0,1 0) 
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Fig.  4.6.  Left:  the  Huygens  principle  at  play  in  the  two  dimensional  case.  Right:  Heap  around  a  circular  cylindrical  obstacle 
for  a  material  with  an  angle  of  repose  of  <5  =  40° 


5.  TECHNOLOGY  TRANSFER.  We  will  pursue  our  collaboration  with  Jenike  &  Johanson,  Inc. 
We  have  had  meetings  about  twice  a  year  and  expect  this  to  continue  in  the  near  future.  The  focus  will 
be  on  the  efficient  resolution  of  the  full  problem  for  granular  flows  (as  described  above),  and  the  fine  tuning 
and  generalization  of  the  powder  consolidation  code.  We  also  intent  to  explore  further  the  possibility  of 
collaborating  with  the  researchers  at  the  U.S.  Army  Engineer  Waterways  Experiment  Station  in  Vicksburg, 
MS,  especially  in  terms  of  the  possible  use  the  present  continuum  approach  to  “plowing  type  problems”. 
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